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ABSTRACT 

We study the internal circulation within the cocoon carved out by the relativistic 
jet emanating from an AGN within the ISM of its host galaxy. Firstly, we develop 
a model for the origin of the internal flow, noticing that a significant increase of 
large scale velocity circulation within the cocoon arises as significant gradients in the 
density and entropy are created near the hot spot (a consequence of Crocco's vorticity 
generation theorem). We find simple and accurate approximate solutions for the large 
scale flow, showing that a backflow towards the few inner parsec region develops. We 
solve the appropriate fluid dynamic equations, and we use these solutions to predict 
the mass inflow rates towards the central regions. 

We then perform a series of 2D simulations of the propagation of jets using FLASH 
2.5, in order to validate the predictions of our model. In these simulations, we vary 
the mechanical input power between 10 43 and 10 45 ergs-sec -1 , and assume a NFW 
density profile for the dark matter halo, within which a isothermal diffuse ISM is 
embedded. The backflows which arise supply the central AGN region with very low 
angular momentum gas, at average rates of the order of 0.1 — 0.8 M Q yr. _1 , the exact 
value seen to be strongly dependent on the central ISM density (for fixed input jet 
power). The time scales of these inflows are apparently weakly dependent on the 
jet/ISM parameters, and are of the order of 3 — 5 x 10 7 yrs. 

We then argue that these backflows could (at least partially) feed the AGN, and 
provide a self-regulatory mechanism of AGN activity, that is not directly controlled 
by, but instead controls, the star formation rate within the central circumnuclear disk. 

Key words: Galaxies: active - galaxies: jets - galaxies: nuclei - methods: numerical. 



1 INTRODUCTION 

There is significant evidence for a close connection between 
the presence of compact objects (hereafter Black Holes, 
BHs) in the dense central regions of most galaxies and the 
AGN phenom enon. Moreover, sin ce the original theoretical 
suggestion bv lSilk fc Reesi (|l99St ). significant observational 
evidence accumulated concerning the impact of nuclear ac- 
tivity o n the global stellar evo l ution within the host galaxy 
(see e.g.lSriiawinski et alj|2006l : lYi et al.|[2007l : iKavirai et all 
12001 120081 . for some recent work). In addition to its effect 
on global star formation, the presence of an AGN seems 
to also be connected with circumn uclear starbursts on small 
(from -parsec to kilo-parsec) scales dStorchi-Bergmann et~al 



200ll : iGonzalez Delgadq fc Cid Fernandes! 1200a ISarzi et a! 
20071 : [TJavies et al]|2007l ). However, despite all this observa- 



tional evidence, the connection between AGN activity, neg- 
ative/positive feedback on star formation, and local star- 
bursts is not well understood. 

If the main physical agent for this connection is the inter- 



action between the jet and the host galaxy's interstellar 
medium (hereafter ISM), then the modeling of the propa- 
gation of AGN relativistic jets into the ISM should have 
as a primary target the understanding of the jet inter- 
action with existing stellar populations, and the mecha- 
nisms which feed the parsec-scale accretion disc around 
the c e ntral BHs. Recent simulations dSutherland fc Bicknelll 
120071 : lAntonuccio-Delogu fc Silk! l2008h have begun to self- 
consistently model the feedback of a relativistic jet on its 
host galaxy, and the global conse quences for e.g. blue -to-red 
cloud migration and downsizing (|Tortora et alj|2009h . 
In this work, we show that the propagation of the AGN jet 
into the ISM not only has an impact on large-scale star for- 
mation within the host galaxy, but also on the feeding of 
the parsec-scale accretion disc around the BH in the very 
central region of the AGN. We show that the generation 
of steep density and entropy gradients near the recollima- 
tion shock (hereafter RS) and the hot spot (hereafter HS) 
induce a backflow within the cocoon, which in turn bends 
back towards the central region, thus providing a secondary 
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accretion flow which is capable of feeding the central BH at 
rates of the order of a few M© yr. _1 , over timescales of a few 
times 10 7 yrs. Thus, this backfiow could easily provide the 
central BH with the "fuel" needed to support the accretion 
disc and the production of the jet, until the RS is destroyed 
and the backfiow disappears. 

This paper is organised as follows. In section 2 we develop 
a model for the origin of the backfiow, by applying an exact 
theorem from fluid dynamics which connects the circulation 
to the entropy gradients (Crocco's theorem). In section 3 we 
describe the results of a series of simulations aimed at vali- 
dating the model developed in section 2, and we also discuss 
the implications for the feeding of the central BH. We draw 
our conclusions in section 4. 

In the following, we express lengths in units of of kilopar- 
secs, but generally, we adopt cgs units, unless otherwise ex- 
plicitly stated. The underlying cosmological paramete rs are 
taken from the 3-year WMAP best fit ACDM model (|Beanl 
l2006h : H = 74 ± 3 Km/sec/Mpc, Q m = 0.234 ± 0.035, 
n b h 2 = 0.0223 ± 0.0008. The unit of time is taken to be the 
Hubble time for this cosmology, i.e.: to — 1.365 x 10 10 years. 
G, kt and m p denote the gravitational constant, Boltzmann 
constant and proton mass, respectively. 



has grown large enough, the jet is pressure-confined and 
the usual scenario applies: the recollimation shock now is 
nearer to the injec ti on p oint, as originally predicted by 
iKaiser fc Alexander! (1 19971 '). The lateral flow now reaches 
the HS and is simply reflectd back, generating a backfiow 
which does not propagate along the bow shock, but is flow- 
ing along the z-axis, tow ards the origin. This has a lso been 
seen in a previous work jKaiser fc Alexanderll 19991 ), includ- 
ing the simulations we present here (Figures [S] and O- 
During both phases, the cocoon is occupied by a very low 
density, high temperature gas. In this region the gas is in a 
turbulent state, with very little or no systematic motions. 
Obviously, the shocked gas within the hot spot has a higher 
entropy than that in the cocoon. Under these conditions, 
Crocco's theorem states that a circulation arise s within a 
compressible fluid, even if in laminar motion fsee lCarJl200ll . 
for a more recent discussion). In a planar geometry one ob- 
tains: 



where: n is the normal to the direction perpendicular to the 
streamline, uj =| uj \ is the modulus of the circulation, and 
s and ho — h + u 2 /2 are the specific entropy and stagnation 
enthalpy, respectively. 



2 FORMATION OF THE BACKFLOW 
2.1 Model 

The model we propose here for the origin and evolution of 
the c i rculat ion within the cocoon extends the models by 
iFallei (|l99lh and IKaiser fc Alexander! (|l997T ). These papers 
show that a recollimation shock (hereafter RS) forms at some 
distance along the path of the jet, when the latter is con- 
fined by the cocoon's pressure. The post-shocked gas then 
accumulates into a terminal region confined between the RS 
and the outer surface of the bow shock(see Figure [TJ, with 
a density n^s ~ 7nj, as appropriate to shocked gas arising 
from a relativistic equation of state. We generically call this 
region the hot spot (HS). 

In the original model by Fall, the only systematic flow 
which is present is that within the jet: the cocoon is assumed 
to be filled with turbulent gas, and the only systematic mo- 
tion within it is assumed to be its expansion within the 
host galaxy's ISM. Moreover, the jet is assumed to be com- 
pletely pressure-con fined by the cocoon. Und e r thes e con- 
ditions, iFaiil (|l99lh and IKaiser fc Alexander] j 19971 ) show 
that a recollimation shock near the injection point of the jet 
develops. 

However, during the initial phases, when the pressure within 
the cocoon is not yet sufficient to confine the jet, the lat- 
ter can expand into the cocoon and acquires a finite open- 
ing angle (see Figure [TJ). We study in detail this phase in 
this paper. If there is an initial transverse velocity gradient 
(dvz/dr ^ 0, the central parts of the jet are supersonic and 
form a shock, while the flow in the outermost part is sub- 
sonic. The shocked gas will accumulate in a terminal region 
called the Hot Spot (hereafter HS). The lateral, subsonic flow 
will eventually also reach the HS and, as it crosses it, will 
acquire an angular momentum and will be deviated back, 
flowing along the inner boundary of the bow shock. 
During a later phase, when the pressure in the cocoon 



2.1.1 Flow near the Hot Spot 

We will first consider the change in circulation induced by 
the presence of the reconfinement shock. We will approxi- 
mate the recollimation shock surface as planar, thus the en- 
tropy is constant across the streamline, and the flow will not 
gain any circulation. However, the gas flowing near and out- 
side the jet will also eventually reach the HS: the boundary 
layer separating this region from the turbulent region will be 
a curved surface, joining the recollimation shock to the bow 
shock (see Fig. [TJ. We then assume that this off-axis flow 
motion is predominantly directed along the z— direction, 
thus: v = v u \ z i. Hereafter, subscripts u and d will denote 
quantities computed in the upstream and downstream re- 
gions, respectively, where by the latter we mean the HS, as 
shown in Figure [TJ In order to compute u>h s , the circulation 
in the hot spot, we will make use of the ex pression f or the 
circulation near a curved shock, derived by IShapirol (|l958l . 
eq. 8): 

2v zW {Ml -I) cos <\> 
Uh3 r s ( 7 + l)M2[2 + ( 7 -l)M2] U 

where: M n — (v z /c s ) cosO is the Mach number in the up- 
stream flow region. As we show in Appendix IB H one can 
obtain an exact solution in 2D spherical coordinates for the 
velocity field after the shocked layer: 

, > _ 4v z \u {Ml - 1) | cos | 
V ^ Ts) - ( 7 + l)Atf2[2 + (7-l)M»] W 

where <j> is a polar angle (see Fig [TJ and, in the coordi- 
nates system shown in Fig. [TJ we have: v z = — v$ cos <j>, v r = 
v,f, sin <j>. We see that the longitudinal z-component of veloc- 
ity after the shock acquires a negative component, and the 
polar component in the downstream region acquires a posi- 
tive value: thus, a backfiow develops. 
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Figure 1. Schematic model of the origin of circulation near the hot spot and in the equatorial plane of the cocoon. We assume that the 
hotspot (HS) is bounded by a lateral spherical section having curvature radius r s . 




Figure 2. The excursion of | v d \ z /v u ^ z |, as measured in eq.[3] for 
different values of the upstream Mach number. The flows within 
the cocoon usually are transonic, with ~ 1 — 2. 



Note that the magnitude of the velocity variation only de- 
pends on 7 and M n . As we can see from Figure [2l this vari- 
ation is not very large, for values of M n ~ 1 — 2, typical 
of the transonic flows wh i ch de velops within the cocoons 
lAntonuccio-Delogu fc Silkl (|2008l . hereafter paper I). 



2.1.2 The flow along the bow shock 

After the HS, the backflow enters into the inner part of the 
bow shock (Figure 0, and the magnitude of the rotation is 
determined by the conservation law for vorticity in 2D (see 
Appendix IA1 for a proof): 



const 



(4) 



Thus, in order to determine | w , we have to determine 
the density inside the bow shock. We notice that this layer 
can not be uniform, because the density, perpendicular to 
the jet, must necessarily be higher than that in the inter- 
mediate regions. In fact, for symmetry reasons, the velocity 
v z m near the meridional plane of the bow shock, thus 
by conservation of the momentum flux, we immediately see 
that the density must increase. 

We will obtain an estimate of the density within this layer 
by making three assumptions: a) the bounding surface and 
its inner part are approximated as oblate spheroids; b) all 
the mass swept out since the formation of the cocoon ends 
in the bow shock; c) the bow shock layer has a constant ma- 
jor axis width Aa. Then the volume of this spheroidal bow 
shock layer is given by: 



AV 



47T 



7? Aa [3aj + 3a;c5a + Aa 



(5) 



where R = bi/ai and at are the aspect ratio and the length 
of the inner semi-major axis, respectively. This gives 



M{a,i) — Anpoa ai 



R 



s/1 -R 2 



arctan(Vl - R 2 ) (6) 



Then, the average density within the spheroidal bow shock 
layer of thickness Aa can be obtained by dividing M(cn) by 
the volume AV: 



p bs = 3po 



A 



3A + 1 



(7) 
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Figure 3. The scaled density within the bow shock as a function 
of a/Aa, the semi-major axis measured in units of the width 
of the spheroidal shell containing the shock. We plot the ratio 
p/(3poa 2 /(Aa) 2 ), f° r * wo diffenmt, values of the aspect ratio R. 



where we have defined A = a;/Aa, and we have dropped the 
subscript i from the semi-major axis. Note that the time- 
dependance of the density enters through the semi-major 
axis (a = a(t)). 

We notice that the evolution of the density in the bow shock 
depends only on the factor A: in Figure [3] the average den- 
sity has a peak at around a fts 0.85Aa, and then decreases 
relatively slowly. 

Knowing the density, we can compute the vorticity of 
the flow after it has been deflected by the HS. We 
will make use of oblate spheroidal coordinates in 2D, 
which define an orthogonal curvilinear coordinate system 
(|Abramowitz fc Steiainll 1972ft : 



x — a cosh £ cos ip 
y — a sinh £ sin ijj 



(8) 
(9) 



As is well known, the scale factors for the (£, ip) coordinates 
are identical: = — a (sinh 2 £ + sin 2 ip) 1 ^ 2 . Curves of 
constant £ are oblate spheroids, with axis ratio b/a = tanh£. 
The eccentricity of these spheroids is then given by: e = 
sj\ - (b/a) 2 = l/cosh£. 

In these coordinates the vorticity u> is given by: 



u> = u, = 



he h 



dip 



(10) 



and the stationary continuity equation can be expressed as: 



V ■ (pv) = 



hch 



dQv^pve) d(h a pv^) 



dtp 



= 



(11) 



Figure 4. The angular factor in eg. 1121 for different aspect ratios 
of the bow shock. 



In Appendix lB2l we use these equations and eq.|4]to deduce 
the tangential velocity v$: 

v l= ■ u2? a) - 2 , [tt + G(0 + 2cosh£-sin 2 ^] (12) 
smh £ + sm ip 

where: B(a) — aco(a)povo, a is a factor depending on the 
initial condition at the injection point, chosen in such a way 
to guarantee that the term inside the parentheses is positive. 
Finally, we have defined: 

G(£) = £Cosh(3£)-~co S h£ 

In the oblate spheroidal coordinate system (£, ip) adopted to 
derive this solution, we have: tpo ^ ip ^ tt/2, and the upper 
limit corresponds to the meridional plane on ellipsoidal sur- 
faces £ = const. We show in Figure [3] the angular variation 
of w.0 , for some values of the eccentricity. As we see, the max- 
imum variations are always around unity, along the stream- 
lines, and the velocity tends to decrease or to stay almost 
constant, as we approach the polar regions. Correspondingly, 
the density must increase or stay almost the same, by flow 
conservation. The variation of the velocity is not very large, 
but depends also on the average density inside the layer. As 
this holds true for all the flow streamlines contained within 
the cross section of the bow shock, we deduce that a high 
density region will form perpendicularly to the jet, near the 
end of the laterally expanding region, for low values of the 
aspect ratio. As a consequence, we expect a steep increase 
of the vorticity for those streamlines approaching the high 
density region, and the flow will consequently bend towards 
the meridional plane. 

In summary, our model makes a few predictions. First, at 
least during the initial phases of the expansion, a circulation 
arises inside the cocoon, induced by the presence of a high 
density and entropy region (the Hot Spot). The backflow 
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which develops tends to follow the bow shock, and then it 
bends back near the meridional axis, perpendicular to the 
jet. For symmetry reasons, this flow will converge back to- 
wards the jet. 

Although we have derived expressions for the velocity along 
a typical streamline (eqs. [3] and I12|l , these should be taken 
to provide only an order of magnitude estimate of the actual 
situation. The presence of an extended turbulent region in 
the central region of the cocoon will introduce fluctuations, 
and we could think of different factors which would make our 
formulae just a first order approximation. Thus it is worth- 
while to use numerical simulations to verify whether this 
model could provide a realistic description of the backflow, 
and the spatial and temporal extent of its validity. 



3 SIMULATIONS 
3.1 Simulation setup 

The observational evidence shows that there exist correla- 
tions between the main physical properties of central BHs, 
AGN properties and properties of their host galaxy halos. 
Some of these correlations, like the relation Mbh — cr, have 
an intrinsic scatter that is relatively small when compared 
to others (e.g. the Pj ct — Mbh relationship) . Our simulation 
setup takes these correlations into account: we have designed 
a set of 9 simulations where the main parameter we change 
is the velocity dispersion of the host DM halo (<x„), and all 
the other parameters are determined by using known scaling 
relations. We determine the tota l virial mass of the h ost DM 
halo using the relation found bv lLintott et al. I (|2006t ). using 
a sample from the SDSS DR3: 

M„ = 2.57 x 10 12 c4 9 9±0 - 15 [km/sec] (13) 

where 0200 is the DM central velocity dispersion in units 
of 200 km sec -1 . M„ is needed to compute the concentration 
para meter of the central NFW halo model, as given by the 
fit of iBullock et~"aT1 (|200ll ): 

(M \-° 13 

Cv = 9 1 w ) ' M * = L5 x wVih ~ lM & ( 14 ) 

We determine the BH mas s Mbh using the relation given by 
iFerrarese fe Merrittl ^OOOD : 

M bh = (1.2 ± 0.2) x lO 8 a^ 7±o ' 3 M (15) 

Finally, we compute t he jet's mechan ical power Pj e t 
using the relation found bv lLiu et alj (|2006h (their eq. 9): 

log(Pj) = -0.22 + 0.59 log (j^J +40.48 (16) 

where Pj is the jet's power in CGS units, and we have 
adopted the value A = Lj, i/L e dd = 0.1 in their equation. 
We then embed an isothermal ISM within this NF W halo, 
follow ing the prescription given in Appendix C of iHesterl 
I2006T ): this halo is specified by the values of the central 
ISM density p C Q and by the virial distance r v . We choose 
the total mass of the hot, diffuse component to be equal 
to the average ratio Ob/fidm appropriate to our chosen 
cosmological model, so that: p;, = 0.174pdm- From this we 
determine the central density p c o and scale radius of the 
diffuse ISM. 




1.0 1.5 2.0 2.5 
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« (Mpc h"') 

Figure 5. Evolution of the velocity field for run slOOav. We show 
the significant transition of the flow near the central region, with 
a radially compressive flow at t rj 13.1 X 10 6 yrs. (top) to a 
more laminar, longitudinal flow at t fa 33.9 X 10 6 yrs. (bottom). 
Flow lines are superimposed on entropy contours. Large devia- 
tions from laminar, null vorticity flows are produced near discon- 
tinuities in entropy, as near the downstream jet's shock in the 
bottom figure. 

We have performed three reference simulations, starting 
with the central values: a v — 100,200 and 300 km sec -1 , 
and with all the other parameters determined according to 
the scaling relations just quoted. In addition, for each of 
these simulations, we have performed two more simulations 
where we have varied M„ (and consequently p C Q and r„ ) 
within the ±lcr bounds, in order to explore a slightly 
larger parameter space. The central ISM density is one 
of the parameters which mostly affect the temporal evo- 
lution of the jet, as recently shown in the simulations by 
ISutherland fc Bicknelll (|2007f ). 
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Table 1. Input parameters of the simulation runs. From left to right, columns 
are as follows: run's label, halo velocity dispersion, halo total virial mass, 
central ISM density, BH mass, jet's input kinetic power. 



run a v M„ p c M bh P k 

(km sec" 1 ) (Mq) (cm" 3 ) (Mq) (ergs- sec" 1 ) 



slOOav 


100 


3.23 


10 11 


2.20 


8.92 x 10 6 


2.29 x 10 44 


slOOpl 




6.44 


10 11 


7.74 






slOOml 




1.62 


10 11 


0.746 






s200av 


200 


2.57 


10 12 


2.68 


1.2 x 10 s 


1.06 x 10 45 


s200pl 




5.68 


10 12 


13.26 






s200ml 




1.16 


10 12 


0.893 






s300av 


300 


8.62 


10 12 


3.04 


5.49 x 10 8 


2.61 x 10 45 


s300pl 




2.03 


10 13 


11.69 






s300ml 




3.66 


10 12 


0.98 
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Figure 6. Same as Fig. 1 for run s200av. Top plot: t Kb 13.1 X 10 6 
yrs., bottom: t pa 44.3 X 10 6 yrs. 



The main input parameters of the simulations are sum- 
marized in Table 13.11 We assume that the jet density is a 
fixed fraction of the central ISM density (pj/pism = 10 -2 ). 
In Table 13.11 we present some quantities related to the 
jet, deduced from the input parameters. It is interesting 
to observe that p c o does not scale monotonicaly with M„ , 
because it also depends on the virial radius r„ . 

As in paper I and in iTortora et all (|2009l . hereafter paper 
II) , all the simulations were performed in a 2D box of size 40 
/i" 1 kpc, with free flow boundary conditions. We have used 
the same setup for refinement and initial mesh distribution 
as in the above cited papers, so we have the same maximum 
resolution: 78.125 pc. The jet is injected with a velocity 
Vj = Vj(y) from the left boundary, with a velocity profile 
given by: 

v * = — uT m ( 17 ) 

cosh{(y-sfe) J j 

where we have defined a dimensionless distance y = r/dj, 
Vi = r c/dj r 'c , being the coordinate of the symmetry ce nter 
of the profile, as we previously did (jTortora et al.ll2009l ). As 
this velocity profile is not truncated, we estimate the size 
of the jet to coincide with the typical scalelength dj. We 
choose dj = 10 fe" 1 pc, and ctj = 10, as previously done by 
iBodo et alj (|l994l ). The peak velocity Vj is determined by 
the total kinetic jet's power: 

dy , , u y (is) 

[cosh{{y -yj) '}] 

We then see that the circulation uj = V x v for this initial 
velocity profile is non-zero. 

The jet is assumed to be non relativistic at injection, as 
the equations we solve ar e also non relativistic. A s it has 
been previously showed bylMarti et al. | (| 1994 1 19951 ) and by 
iKomissarov k, Falld (|l996t ). it is not possible to find an ex- 
act relativistic anologue of a classical jet, but one can find 
configurations which are physically comparable. We can de- 
rive the main parameters of a relativistic jet which has the 
same kinetic power and radiu s as our jet, using eq. (21) from 
ISutherland fc Bicknelll (|2007l ): 

P k , r = 3.9 x l0 40 -^ I ^ 7 r r d? r j 2 ft (l + ^^x) (19) 



where 7 is the adiabatic index of the ISM, £ = pj/pi 3rn , 7T7 = 
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Table 2. Main parameters of the jet, deduced from the input values. From left 
to right, columns are as follows: run's label, (spatial) average input velocity, 
Lorentz factor, relativistic velocity, (non-relativistic) Mach number, enthalpy 
flux, and ratio between enthalpy and kinetic flux. Note that for an ideal gas 
the enthalpy is given by: h = "fp/(-f — 1), thus the enthalpy provides also a 
measure of jet's pressure. 



run 


<v>, 




Pi 






Ph/Pk 




(km-sec ) 








(10 16 ergs- sec 1 ) 




slOOav 


217.4 


6.74 


0.989 


11.55 


22.83 


0.029 


slOOpl 


143.0 


3.65 


0.961 


6.08 


80.1 


0.108 


slOOml 


312.8 


11.55 


0.996 


19.93 


7.72 


0.010 


s200av 


340.7 


13.11 


0.997 


22.65 


27.74 


0.008 


s200pl 


199.3 


5.93 


0.985 


10.12 


137.3 


0.039 


s200ml 


491.0 


22.68 


0.999 


39.25 


9.25 


0.003 


s300av 


439.2 


19.27 


0.998 


33.33 


31.44 


0.004 


s300pl 


281.0 


9.84 


0.994 


16.95 


121.05 


0.014 


s300ml 


64.1 


33.9 


0.999 


58.72 


10.14 


0.001 



(? w/fc)/10 7 , dw = (2d J )/(Wpc), X = PjC 2 /A P] , and Y 3 = 

1 / .J 1 — P | is the jet's Lorentz factor. Given as input Pk, r 

and dio = 2, we compute the Lorentz factor (see Table l3~Tj) . 
We also compute the ( non relativistic) Mach number , using 
eqs. (17) and (18) from ISutherland fe Bicknelil (|2007T ): 

M tH . = ^ I (l+^ X )(r?-l) (20) 

We note from Table 13.11 that the jet is hypersonic at injec- 
tion, and that the mass flux is very small. As we will see, 
the mass flow from the backflow is significantly larger: most 
of the flux contributing to the latter comes from ISM gas 
dragged by the jet itself. 

The total energy flux injected by the jet is given by: Pj = 
Pk+Ph, where Pk was given above (eg. ITS)) , and the enthalpy 
flux is given by: 

P h = ird}h{Vj) (21) 

where h is the specific enthalpy, which in a fully ionised jet 
is given by: 

h = p (22) 

7-1 

We see from the last column of Table 2 that the enthalpy 
flux is comparable and even much larger than the kinetic 
flux. Thus the entrained ISM gas gets a significant thermal 
energy input from the jet, in addition to the kinetic drag. 

3.2 Simulation results 

The simulations confirm the general picture of our model. 
Two different circulation regimes descr ibe the flow within 
the cocoon . As in itially demonstrated bv lFalld |l99ll ) and by 
lAlexandeil l|2006l ). the jet forms a recollimation shock along 
its path, when the average pressure p c inside the cocoon is 
comparable to the jet's lateral pressure. The jet then is con- 
fined towards the axis. Its injection flow is supersonic, and at 
its terminal point, near the hotspot, the jet moves subsoni- 
cally: thus, a shock develops. However, the simulations show 
that the position and extent of this shock changes with time. 
As the character of the backflow is critically dependent on 
the shape and position of this recollimation shock, we will 



analyze here in some detail this issue. As we can see from the 
upper parts of Figs. [5] and [6] initially this shock is perpendic- 
ular to the flow and has a transversal extent comparable to 
the cocoon's extent: also, the post-shocked gas flows along 
the axis, thus the gas does not gain any circulation. How- 
ever, the gas flowing off-axis eventually encounters the den- 
sity discontinuity near the HS, and starts to gain rotational 
motions. This then generates a backflow, which transports 
gas in the direction opposite to that of the inflowing jet. This 
backflow then changes direction again when it reaches the 
meridional plane, near the vertical axis of the configuration. 
Note that all the four boundaries here are not reflecting: the 
fluid can freely leave the box, thus the change of flow direc- 
tion is not a result of particular boundary conditions. The 
resulting flow could then result in compression of a gaseous 
disk perpendicular to the jet and lying in the meridional 
plane. 

The reconfinement shock is however unstable: as soon as 
the cocoon expands, transonic turbulence develops which 
destroys the shock. The jet is now confined by the turbu- 
lent pressure of the cocoon, but the extent of the recollima- 
tion shock is much smaller, because the cocoon's pressure 
p c steadily increases, as the jet continues to input thermal 
energy within the cocoon, whose expansion is slowed down 
by the ISM pressure. Now the pattern of the flow changes 
drastically: as we can see from the bottom parts of Figs. 
and[S] the flow is mostly directed horizontally towards the 
accretion disc. The recollimation shock has now a much 
smaller transverse extension, and is loca ted very ne ar to 
the injection point, a s orig inally shown by|Fallc (199l]) and 
iKaiser fe Alexander! (| 19971 ). This general evolutionary pat- 
tern of the flow is reproduced in all the simulations we have 
performed, although with some variations in the temporal 
evolution and some large transient fluctuations (see Fig. [8}. 
In Figure [7] we present a schematic summary of the result- 
ing backflow around the circumnuclear disc, during the two 
main phases just mentioned. During Phase 1 the backflow 
mainly comes along the meridional plane, thus its angular 
momentum is almost zero, while during Phase 2 the back- 
flows has a finite angular momentum component out of the 
plane. The part of these backfiows falling outside the central 
nuclear disc will collide and lose most of its angular momen- 
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Backflow 



Phase 1 



Backflow 




Phase 2 



Figure 7. Schematic view of the backflows during the two phases. We show a model of a realistic situation, when a two-sided jet/cocoon 
develops, and there are backflows from opposite directions. We assume that the circumnuclear disc has a finite radial extent, charcatcrizcd 
by a finite radial scalelength l^. 



turn, flowing then radially towards the dies, as in Phase 1. 
The final fate of the backflow during Phase 2, when it falls 
on opposite sides of the disc, will depend on whether the lat- 
ter is mostly stellar or gaseous: we comment more on this in 
section l3.3.1l During both phases, most of the backflow will 
then result into a net flow of almost zero angular momentum 
gas towards the central regions of the AGN. In the next sec- 
tion we compute the actual fluxes around this circumnuclear 
region. 

3.3 Inflow around circumnuclear regions 

In so much as there exists a shock along the jet, a fraction 
of the backflow towards the z~0 plane is present. This is 
true both during the initial phase, when the shock is rather 
extended, and also later, when the flow is reflected by the 
terminal shock. This backflow mostly entails gas supplied by 
the lateral, subsonic jet flow. The general circulation induced 
by the backflow contributes a systematic flow of low angular 
momentum gas in a plane orthogonal to the jet, near the 
z=0 plane. It is then interesting to compute the mass flux 
contributed by this backflow to the very central region, near 
the central Black Hole. 

If an accretion disk is present around the BH, it would lie 
in a plane perpendicular to the jet, and its structure will 
probably be affected by the backflow. In Fig. [9] we show 
the average mass flow within an annular region having an 
internal radius ri n — 1.5 pc, a radius r c — 20 pc and height 
h c — 2 pc, for the three main runs. The mass flux across 
the lateral surface of the disc will be given by: 

M r = 2irr c / dzpfv r (23) 
Jo 



assuming axial symmetry. There will also be a flux across 
the upper surface of the disc: 

rr c +d c 

M~ =2ir drrpfv z (24) 

We show in Figure [10] the evolution of the total mass flux 
M r + M z . It is interesting to observe that the accretion rates 
contributed by the backflow all seem to vary on similar time 
scales, of the order of 3.5 — 4 x 10 7 yrs., during which the 
accretion rate can reach peak values in excess of 1 Mq yr" 1 , 
i.e. the typical accretion rates required to support activity 
in QSO and Seyferts (see e.g. Ijogeel \20o4 . Table 1). This 
time-scale seems to be almost independent of the ratio be- 
tween the bulge and BH masses, but we see large variations 
induced by variations in the local ISM density (Fig. I10[) : 
systems where the local ISM density is larger seem also to 
have larger accretion rates, with peak values exceeding 4-5 
Mq yr~ J . This suggests that the accretion rate can be very 
sensitive to the occurrence of recent episodes of wet mergers, 
during which the gas density can easily vary by ±lcr w.r.t. 
the average values. 

It is also interesting to observe that the accretion is highly 
intermittent, as can be appreciated by looking at Fig. [8] 
The beginning of the intermittency coincides with the dis- 
appearance of the reconfinement shock, and marks the onset 
of fully developed, transonic turbulence within the cocoon. 
If part or all of this backflow feeds the central BH, the in- 
termittency could then translate into an intermittency of 
the feeding mechanism, and result in a random orientation 
of the jet axis with the structure of the surrounding galaxy 
(jKinnev et al.ll2000h . 

Finally, we have also plotted in Fig. [l0]the predictions about 
the mass flow rates from the model that we developed in the 
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Table 3. Mass inflow rates. We show the input mass flow rate, (M); n (column 
2), and the (averaged in time) flow entrained within the backflow around the 
circumnuclear region (column 3). Flows are in units of Mqjt. - 1 . 
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Figure 8. Mass flow in the circumnuclear region for run 
sim200av. We show the averaged flux (continuous line) and the 
instantaneous flux (dotted line), both computed by adding the 
contributions from eqs. 1231 1241 as described in section [3.31 . The 
transient fluctuations are caused by the intermittent behaviour of 
the flux inside the cocoon, during the evolved phase of the flow. 



previous sections. The model seems to adequately describe 
the results from the simulations. The quantitative agreement 
is quite good, at least during the phase dominated by the 
regular backflow, i.e. the only phase which we have mod- 
elled. We have adopted the temporal dependence u sing the 
expansion law given bv lKaiser fc Alexander! (1 19971 ). for the 
input values of Pk and rii srn as given in Table 1 (we con- 
sider only the fiducial models). We notice that the largest 
deviations are for smaller values of a v , when the expansion 
is slower and the deviations from the self-similar model by 
Kaiser & Alexander are larger. 

In Table 3 we compare the input mass flow with the flow 
entrained by the jet: the latter is about two orders of mag- 
nitude larger than the first. The combined action of kinetic 
and enthalpy pressure exerted by the jet on the ISM is able 
to create a circulation which drags a significant amount of 



Figure 9. The average mass flow rate around a circumnuclear 
region, with a diameter d = 20 pc centered around the circumnu- 
clear disk. The flow is computed as described in Figure [8] 

We show only the results for runs sSOOav (continuous line), 
s200av (dotted line), and sSOOav (dashed). Time scales are 
almost constant, while the maximum of the curve seems to 
depend on P^. 

mass flow towards the central regions. Note that this how- 
ever takes place on a time scale of the theorder of few times 
10 7 yrs. 



3.3.1 Feeding the central BH 

We now briefly consider the final fate of the inflowing gas. 
As we have seen, this gas has very low angular momentum, 
at variance with most of the gas in the outer regions of 
the circumnuclear disc. However, the final fate of this gas 
will primarily depend on the structure of the circumnuclear 
disc around the central BH, and in particular on whether 
this disk has a pressure which can signficantly resist the 
inflow . Since the seminal paper by IShlosman fc Beeelmanl 
1 19891 ). it has been suggested that on scales larger than 
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Figure 10. Mass inflow rates around a central disk region, com- 
puted as described in Figurc|8] The continous lines are for simula- 
tions s(100,200,300)av, dotted lines for s(100,200,300)pl, dashed 
lines for s(100,200,300)ml. The blue dashed lines are the fits from 
the model given in the text. The red vertical lines mark the epoch 
when the reconfinement shock is destroyed and the jet starts to 
propagate freely, being only confined by the cocoon's ram pres- 



10-50 pc, the circumnuclear discs around AGNs are prone 
to be unstable to gr a vitational fragmentation (iBegelmanl 
1994 lGoodmanll2003l; IWada fc Normanl |2002|; IWadall2004 



Tan fc Blac kman 2005; Escala 2007; Es cala fc Larsonll200ct ) 



Recently, it has been suggested that star formation within 
this disc can provide a self-regulatory mechanism which has 
a two-fo ld action: heat i ng the disc, thus halting star for- 
mation (| Collin fc Zahnl Il999al lbl) , a nd feeding the central 
BH through massive stellar winds (|Thompson et all 120051 ; 
iNavakshin et all 120091 ). More recent observations with res- 
olution on parsec scales are now confirming the presence 
of dense star clusters i n the central parts of bulges which 
are also hosting AGNs |Davies et al.ll2007l: ISeth et al.ll200Sl ; 



iGonzalez Delgado et alj 120081 ; iGraham fc Spitlerl 12009? ) . 

thus providing some evidence for the existence of gas and 
especially stars within these circumnuclear discs, from a few 
to a few hundreds of parsecs. 

Thus, we expect that the inflowing gas provided by the large- 
scale backfiow will mostly "see" a stellar disc, and will inter- 
act with the stars by creating local Bondi-like accretion flows 
around them. The cross sections for Bondi flows are very 
small (of the order of few stellar radii), and we can then rea- 
sonably predict that most of the inflow will be able to reach 
the inner gaseous, high-density (n c > 10 8 cm -3 ) disc, thus 
providing significant feeding of the latter. In a forthcoming 
paper (Antonuccio-Delogu & Silk 2010, in preparation), we 
will analyse this interaction in detail to understand quan- 
titatively how much the inflow can contribute to actually 
feeding the AGN. 



4 CONCLUSIONS 

The development of a backfiow within the cocoon that 
results from the interaction of a relativistic AGN jet 
with the surrounding ISM has been previously studied by 
iKaiser fc Alexander! l fl999h . Our main focus in this work has 
been to show that a significant contribution to the develop- 
ment of circulation within the cocoon comes from the pres- 
ence of naturally developing density and entropy gradients. 
Our calculation can have direct implications for the struc- 
ture of AGNs' and of their circumnuclear discs. One of the 
best studied cases is t he ionised circum nuclear disc in the 
central region of M87 (|Ford et alj|l994 ). HST spectroscopic 
observations have s hown show that this disc is shock-excited 
|Dopita et al.lll997h . and its emission lines are probably orig- 
inating from accretion of ga s external to the di sc (see how- 
ever the critical remarks by ISabra et al.l (|2003f 0. It is then 
possible that the backfiow we have studied in this work could 
contribute to the accretion shock. We plan to present a de- 
tailed modelling of the contribution of this flow to the pb- 
served spectrum in a future work. 

One of the first questions concerns the stability of this flow, 
and a useful hint towards an answer comes from the simula- 
tions we have performed. Although a linear stability analysis 
would also be possible, we see that the flow is not heavily 
affected by fluctuations induced by turbulence within the 
cocoon. However, when the recollimation shock is destroyed 
by instabilities within the cocoon, the main source of cir- 
culation, which arises from discontinuities within the HS, is 
destroyed, and the backfiow changes abruptly, becoming al- 
most totally longitudinal. The shape of the jet now changes: 
the recollimation shock forms again, but much neare r the 
injection point, as o r iginal ly predicted by iFalld (|l99lf ) and 
IKaiser fc Alexander! (1 19971 ). and a new shock forms imme- 
diately before the hot spot, from which now a backfiow de- 
velops which flows almost parallel to the confined jet. Also 
this backfiow contributes very low angular momentum gas 
to the central regions of the AGN, with average mass inflow 
rates of the order of a few M© yr. , on time-scales of the 
order of a few times 10 7 yrs, as required in order to feed an 
AGN. 

In summary, we can distingusih two main phases in the de- 
velopment of flow within the cocoon: 

• Phase 1 : A large scale backfiow develops, fed by some 
gas originating in the lateral regions of the Hot Spot, and 
flowing along the inner region of the contact discontinuity, 
up to the meridional plane of the cocoon, and eventually all 
the way down towards the AGN; 

• Phase 2: When the jet becomes pressure confined, gas 
originating from the HS flows back longitudinally towards 
the AGN. 

During both phases the backfiow provides (almost) zero an- 
gular momentum gas to the central region of the AGN. 
We have seen that the backfiow originates when the subsonic 
flow of the outer parts of the jet hits some shock, and then 
gains angular momentum. Thus, this flow is not arising from 
some turbulent entrainment, as is seen when part of the 3D 
turbulent energy is dissipated through shear and engenders 
a regular motion of entrained gas. This is the reason why 
this backfiow is also manifest in 2D simulations, as those we 
have studied in this work. 
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Competing mass inflow mechanisms, in particular low angu- 
lar momentum gas driven inwa rds by a nuclear s piral disk, 
can power nuclear starbursts (| Fat hi et al.l l200Sl ) but give 
substantially sm aller mass inflow rates within the central 
few parsecs (e.g. iDavies et alj|2009f ). 

It is interesting to observe that the inflow rates appear to 
depend strongly on the local ISM density: more massive 
and denser galaxies host stronger backflows, because there 
is more matter which is dragged in by the jet. The average 
mass inflow rates in the circumnuclear regions are always 
of the the order of 1CP 1 M sun yr. _1 , although they do not 
scale only with the jet's input power, but seem to be strongly 
affected by the ISM density, as can be seen from Table 3. 
The kinetic pressure of more powerful jets is typically much 
larger than the ram pressure from the gas: only within the 
turbulent region of the cocoon, where the backflow is absent, 
does one reach approximate equipartition (paper I). 
This backflow could thus provide the main source for feed- 
ing the central AGN with very low angular momentum gas, 
which would eventually reach the BH and feed the jet. We 
see that, when this backflow is present, the jet power in- 
creases with time, and the cocoon expands more rapidly. 
Ultimately, this would have the effect of accelerating the 
expansion of the cocoon and destroying the backflow, thus 
halting the feeding of the AGN. From Figure 9 we see that 
the mass inflow rate initially increases, thus providing more 
fuel to the AGN, but then decreases, as the backflow now 
originates from gas flowing near the hotspot, whose energy 
decreases as the cocoon expands. One can presume that a 
fraction of the inflowing gas will be able to reach the AGN 
and feed it, possibly increasing Pk. This however depends 
on the fate of the inflowing gas when it meets the inner ac- 
cretion disc and the stellar disc, i.e. on phenomena taking 
place on parsec scales, much smaller than those resolved by 
the simulations we have presented in this paper. 
We have thus described a self-regulation mechanism, which 
does not directly invoke any starburst in the circumnuclear 
disk to feed the AGN. The role of starbursts in AGN feeding 
is that of supplying low angular momentum gas from AGN 
winds, which can eventually reach the central black hole. In 
this backflow model, the low angular momentum gas is pro- 
vided from the jet itself: however, more detailed modelling 
of its interaction with the accretion disk is needed to under- 
stand quantitatively how much this inflow can contribute to 
the feeding of the AGN. 
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APPENDIX A: A VORTICITY CONSERVATION 
THEOREM FOR 2D MOTION 

In the particular case of a 2D motion, when the circula- 
tion u) is everywhere perpendicular to the velocity field, it 
is possible to show that the quantity | u) \ /p is constant. 
For a general 2D co mpressible motion the evolution of the 
vorticity is given by (|Colelll962l , eq. 1.46): 



dui „ 
— = -uA7 ■ 

dt 



(Al) 



where, as usual, the total time derivative denotes the con- 
vective derivative. 




Figure Bl. Coordinate setup used to solve the circulation equa- 
tion near the HS region. 



From the equation of continuity we also have: 
dp 



dt 



-pV • v, 



thus: d\np/dt = — V ■ v, and, after substituting in eq. IA11 
we eventually obtain: 

duj d In p 
— = u) 

dt dt 

Dot multiplying both sides by lj, we arrive at: 

1 d 1 i |2 d 1 i i d 1 
In w = — m \ ijj \— — In p 



2dt ' 1 dt 
which can be rewritten as 



dt 



H^ 1 



(A2) 



Thus, the quantity | w | / p is advected along the streamlines. 



APPENDIX B: 
EQUATIONS 



SOLUTIONS OF FLOW 



We will here describe a scheme to derive approximate solu- 
tions of the flow equations along the backflow. As described 
in the text, we distinguish three regions in the backflow: a) 
the Hot Spot, b) the bow shock, c) the meridional region. 
We will now shortly describe the solutions in each of these 
three regions. 

Bl Streamlines across the HS 

We make reference to the geometry described in Figure [BTI 
We define the polar angle <f> (not shown in figure) as the angle 
between the r-axis and the position vector £. Introducing 
an auxiliary angle 9, we see from the figure that a — n — 
(tt/2 + 9) = tt/2 - 9, and: <j> = tt + (tt/2 - 6) = 3tt/2 - 9; 
thus, we have: a = <f> — 7T. We will find useful to compute 
the circulation in a reference frame centered at the origin 
of the (assumed) spherical region bounding the sides of the 
HS. In a 2D spherical coordinate system, the out-of-plane 
component of the circulation is then: 



(?t^| d ) 



The continuity condition across the normal gives: 

PuV z \ u cos9 = pdV£\d 



(Bl) 



(B2) 



where the subscripts "u" and "d" denote quantities evalu- 
ated in the up- and downstream regions, respectively. After 
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substituting v^ d from the latter equation in eg. IB1I and as- 
suming that the downstream region is homogeneous, we see 
that only the first term in the right-hand-side of eq. IB1I is 
different from zero. 

We now apply Crocco's theorem, making use of the expres- 
sio n for the circu lation after a curved shock layer, derived 
bv lShapird j 19581 , eq. 8): 



{Ml 



«(<ft) 



£ d£ r s (7 + 1) Ml [2 + ( 7 - 1) Ml] 



(B3) 



We now integrate both terms of this equation across a layer 
of width Ar, starting from the layer at £ = r s , thus obtain- 
ing: 

(r s + Ar)u 0!d (r 3 + Ar) - r s v^,\ d (r a ) = 



(2r s + Ar) Ar- 



1) | cos(<ft) 



(7 + l)M2[2 + ( 7 -l)M2 



(B4) 



Passing to the limit Ar — > 0, after having expanded v^,(r a + 
Ar) in Taylor series, we finally obtain: 



V4,(r s ) + r. 



5=r s 



4v zlu {Ml - l) | cos </> | 
(7+1) Af2[2 + ( 7 -l)M2] 



(B5) 



As it is evident from this equation, we have yet to specify 
the gradient of the angular velocity across the boundary: 
however, we will now show that, by applying the continuty 
equation, the latter condition can be replaced by an initial 
condition on the circular velocity v^ . Until this point we 
have only used the equation for the circulation, but we can 
get an additional constraint using the continuity equation: 

V ' ^ = £ 8i ^PdV^a) + ^qJ {PdV^\d) = (B6) 

Inserting in the latter equation the expression for pdV^ from 
eq, IB2l and using the homogeneity in the downstream HS re- 
gion (dpd/d(j> = 0), we eventually get the following equation 
for the angular velocity: 



1 ■ j. , P ddv <f> n 

7 • Puv z , u smtft + — — — = 

which immediately provides the following solution: 



v 4>\d = —v z \ u (1 + cost 

Pd 



+ v 4> (£></> = ?>") 



(B7) 



(B8) 



Note that cos<ft = — sin#, thus, imposing that eqs. IB7I 
and IB8I describe the same solution, we obtain: 



P„ 
Pd 



Mi 



1) sin(6>) 



7 + 1 M2[2 + ( 7 -l)M2] 



and: 



V4. (£, 4> = t) = -v z \u + r. 

pd 



into the continuity equation 
latter to obtain: 



. 1111) we can integrate the 



h^pv^ - h^pv^ \^ = ^, =p(£) 



(B9) 



where p(£) is an integration constant. Using this latter equa- 
tion, we can the express the density as: 



g(0+p(0 



(BIO) 



where we have defined: q(£) = h^pv^, |v>=v>o- We now make 
use of the law of vorticity conservation (eq. IA2|) to deduce 
that: 

1 d 



(/i^ity,) = c (a)p 



(Bll) 



where the function Co (a) is determined from the initial con- 
ditions, i.e. from the typical vorticity and density at the 
injection point. Substituting the density from eq. IB10I in 
eq. IB 1 1 1 we obtain: 



1 d 



(/i^iy,) = Co (a) 



g(0+P(0 



This can be integrated to give: 

h\v% = h\v% | e = 5o +a 2 [«(£) + {«(£) - v(£ )} sin 2 ip] -a 2 w(£ ) 

(B12) 

where we have made usage of the expression for the scale 
factor h^ t ii> given in section[2] and introduced the definitions: 



-(0,0 = 00(0) [?(£)+ P(01 



and: 



u(£) = 2 J dar(a) sinh 2 cr, «(£) = 2 j dar(a) 

p(£,,ipo), and: vo = v^(£,rpo), we see 



Now, defining: po 
that: 



.1/2 



?(£) = PoVoa (sinh 2 £ + sin 2 ip Q ) ' 

Without loss of generality, we now assume: p(£) = 0. Thus: 
r (£) = c o( a )9(£)i and, subsituting in the expressions above, 
we have: 

1 

6< 

and: 



u(0 



-copovoa [cosh(3£) — 9cosh(£)]f. 



(B13) 



= 2copov acosh(£) l]i 



(B14) 



In deriving these expressions, we have assumed that 
sin 2 ipo ~ 0, i.e. that the injection point of the bow shock is 
very near to the equatorial plane, a reasonable assumption 
if the Hot Spot lies near the tip of the expansion zone. In- 
serting eqs. lB13l and lB14l into the general expression eq, IB12l 
we eventually get the solution: 



2 



acopovo 



sinh £ + sin ip 



B2 Solution along the bow shock. 

To describe the motion in the region of the bow shock, we 
will make use of a more convenient oblate spheroidal coordi- 
nate system, as defined previously in section [2] We will also 
assume that 3> V(, i.e. that the expansion of the cocoon 
is negligible w.r.t the meridional circulation. Inserting this 



1 3 

— cosh(3£) — — cosh £ + 2 cosh £ • sin ip 



(B15) 



where, without loss of generality, we have assumed that: 
h\v% | fj =0. 




ram pressure 



